Noise Effects in Nonlinear Biochemical Signaling 



Neda Bostani/ David A. Kessler,^'0 Nadav M. Shnerb,^'|^ Wouter-Jan Rappel,^'[^ and Herbert Levine^'j^ 

'Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, 
Chinese Academy of Sciences, Beijing 100049, China 
^Department of Physics, Bar-Ilan University, Ramat-Gan 52900 Israel 
^Center for Theoretical Biological Physics, University of California San Diego, La Jolla, CA 92093-0319 USA 

It has been generally recognized that stochasticity can play an important role in the information 
processing accomplished by reaction networks in biological cells. Most treatments of that stochas- 
ticity employ Gaussian noise even though it is a priori obvious that this approximation can violate 
physical constraints, such as the positivity of chemical concentrations. Here, we show that even when 
such nonphysical fluctuations are rare, an exact solution of the Gaussian model shows that the model 
can yield unphysical results. This is done in the context of a simple incoherent-feedforward model 
which exhibits perfect adaptation in the deterministic limit. We show how one can use the natural 
separation of time scales in this model to yield an approximate model, that is analytically solvable, 
including its dynamical response to an environmental change. Alternatively, one can employ a cutoff 
procedure to regularize the Gaussian result. 

1^ PACS numbers: 02.50.Le, 05.65. -(-b, 87.23.Ge, 87.23.Kg 

l> 

I I I. INTRODUCTION 

The role of stochasticity in the functioning of cellular signal transduction networks is a question of great topical 
interest [1]. Unlike typical condensed-matter systems, biological cells must carry out chemical manipulations with 
O small numbers of molecules, an inherently noisy situation. Noise comes in a variety of forms, including fluctuations 
in chemicals to be sensed [5], fluctuations in the binding-unbinding of receptor arrays [3], fluctuations during the 
processing of information [3] , and fluctuations in the implementation of downstream actions [5] . 

In this context, almost all analytic studies of stochastic reaction dynamics utilize a small noise Gaussian approx- 
imation. This assumption emerges naturally, for example, in the van Kampen system-size expression [5] where the 
^ fluctuations in particle number are formally lower order and hence are treated as small and centered around the 
^ mean value set by the deterministic reaction equations. The initial purpose of this paper is to point out that this 
•/"J approach may give highly misleading results especially when some of the downstream reactions are nonlinear. We do 
this by studying a speciflc example, that of an incoherent feedforward module processing data from a small number 
of receptors [7]. Afterwards, we show how an alternate approach for the rapid activation versus slow inhibition limit 
can provide a complementary analytic approach. 

The example we choose to study is a part of the gradient sensing module underlying the chemotactic response of 
Dictyostelium cells [5j . These cells appear to implement a control circuit incorporating a simple incoherent feedforward 
^ I loop topology for adapting out the constant concentration background [51 . This circuit is instantiated by using 
. . a RAS-GEF as a positive signal and RAS-GAP as the complementary inhibitor [TT]. These are both activated upon 
^ the binding of cAMP by the G-protein coupled cAMP receptor and in turn drive the signaling hub protein RAS into 
k>( its active (GTP-bound), respectively inactive (GDP-bound) form. The circuit diagram is illustrated in Fig. [l] In the 
limit where we are far from saturation, the system can be described by the following equations [10; : 

f .,S(0-» 
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FIG. 1: Circuit diagram showing the activation of both A and B by the signal S. A in turn activates E, while B inhibits it. 

Here E is the fraetion of RAS moleciilcs that have been activated and S is the external signal which drives both 
the activator A and inhibitor B. It is trivial to verify that if 5 is a constant signal, the steady-state value of E, 
Eq = Q.6/{a5 + (3^) does not depend on its value. Hence this system in the deterministic limit is a perfectly adapting 
module, exhibiting only a transient response to changes in its input. 

In order to study the effect of noise in the input signal S{t) on this system, it is standard to assume that S{t) is 
the sum of a deterministic signal So{t) plus stationary Gaussian noise. For simplicity, we consider the case where the 
deterministic signal is a constant, Sq, and the noise, r]{t), has zero mean and correlator 



-\t-t'\/r 



(2) 



The advantage of assuming Gaussian noise is that the system is then analytically tractable. The fields A and B can 
be expressed in term of the noise rj as follows: 



A{t) =^So + a f dt' 
B{t) =f5o + /3 f dt' 

J — OO 



-i(t-t') 



-6(t-t') 



Substituting this into the effector equation allows us to find 



E{t)= f dti5i(ti)e"^*'i'*^''^^*=^ 

J —OO 



(3) 



(4) 



with the definitions 
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5o+ f dt' (ae-T(*-*')+/3e-*(*-*'))r/(i') 

J —OG 



(5) 



From this expression, all moments of E can be calculated exactly. For example, let us focus on the expectation values 
of E. The standard expressions for Gaussian processes, for example 



- J{S{t')-s,g)h{t')dt' 



)=e-^ 



^ J dt'dt"h{t')G{t',t")h{t") 



(6) 
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FIG. 2: (color online). Left) Excerpt of a simulation of the Gaussian model with So = 1/3, = 2/9, a = 1, 7 = 2.5, 5 = 0.4, 
/3 = 1.7, showing a large negative fluctuation of E. Note that /3 is chosen to be exactly the critical value at which Eq. (|9| 
is violated and {E) diverges. Also shown are S{t)/So, A{t)/Ao and B{t)/Bo. Right) Larger time series of E{t) showing the 
intermittent nature of the large fluctuations. 



allow us after a tedious calculation to derive the following expression for {E) 

/o V 7 



(7) 



From this, one can immediately recover the aforementioned a'^ — Q deterministic result, {E) — gj"p^ ■ The exact 
expression for Ae{u), given in the appendix, is not particularly informative; the only critical feature is that it decays 
to a constant at large u. The factor cr^$(u) is the exponent is much more significant. Again, we leave the full form 
for the Appendix, and merely give the large u behavior: 



a B 

- + 7 
7 



(8) 



This term represent the diffusive growth of (the integral of) 5^ and has a well-defined form even in the white-noise 
limit for S where r — !■ with cr^T fixed. 

The starting point of our work is the observation that the integral defining {E) fails to converge unless 



5o > crV 



(9) 



ft is easy to show that similar, but more stringent, bounds hold for all moments of i?, which therefore are predicted 
to grow without bound (starting from any initial state) if the noise is too large. The problem arises from the fact that 
the S fluctuations are unbounded below for Gaussian noise and hence can drive A and/or B negative. This gives rise 
to transient periods during which \E\ grows exponentially. This behavior can be directly seen in a simulation of the 
Gaussian noise model, as presented in the left panel of Fig. [2j (In passing, these large excursions are rare events if 
tr is sufficiently small, and hence getting an accurate measure of {E) from the simulations can be difficult. We will 
return to this point in more detail later). As the noise gets large these growth periods lead to ever increasing values 
of \E\ and the anomalous contributions to {E) never saturate. It is important to note that even when Sq/u <C 1, 
so that the negative fluctuations of A and B are rare, nevertheless the Gaussian model can yield unphysical, infinite 
results, if a/7 + /3/(5 is large enough. Thus, our finding depends essentially on the nonlinear coupling of the noisy 
signal to the E field, so that the noise acts multiplicatively on E. 



II. BINOMIAL NOISE 



Clearly, the Gaussian noise approach is in general unacceptable. We must treat the noise in a more realistic fashion 
if we are to have a well-defined model. If the source of the noise is the finite number of receptors [TH |T3] , we are led to 
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consider a model wherein the incoming signal is a random variable reflecting the fraction of bound sensors. Assuming 
that the N receptors are independent and bind a ligand of fixed concentration, cq, the signal can be exactly described 
via the master equation for the probability distribution, P(s, t), for the number s = 0, 1 . . . , iV of occupied receptors, 
where the signal S — s/N] i.e., the fraction of occupied receptors. 

dP(s t) 

g^' ' = -[fcoffS + fco„co(^-s)]P(s,t)+fcoff(s+l)P(s+l,0 

+ KnCo{N - S + l)P{s-l,t) 

We call the model with this discrete noise the binomial noise model, as the equilibrium distribution of s is binomial, 

o , ^ _ N\ (fconCo)-(fcoff)^-- 

For large N , we recover a Gaussian noise process, except in the tails, as can be seen via the following argument. Our 



discrete stochastic process, Eq. [TOj is well approximated for large N by the Ornstein-Uhlenbeck process, described 
by the Fokker-Planck equation 

with So = fconCo/(fconCo + fcoff) , = fcoff + fconCQ and cr| = fcofffconCo/(iV(fconCo + koff)'^) SO that in the limit the 
stochastic function S{t) is Gaussian with variance as- Here, of course, the equilibrium distribution is 

PGMS)c^exp[-{S-Sof /2<jI'\ (12) 

and the steady-state autocorrelation function is 

< S{t)Sit') > - 5*0 = G(<, t') = g-|t-t'|/r (13) 

Nevertheless, for all finite N, the signal S in the binomial process is always non-negative, and no anomalous behavior 
can occur, with E(t) strictly bounded from below by 0. The striking difference in the two models is apparent by 
comparing the Gaussian simulation presented in Fig. [2] to the simulation of the corresponding binomial model in Fig. 
[3l The well-behaved nature of the binomial model for all N is consistent with the fact that our condition for the 
convergence of the first moment in the Gaussian model, Eq. ([9), is always satisfied in the large N limit, given that 
the noise amplitude is small, of order 0{1/N). The large- limit, however, is only an asymptotic approximation, 
since for any finite value of the noise amplitude cr| of the Gaussian model, sufficiently high moments of E do indeed 
diverge. 

An additional perspective on the difference between the binomial and its parallel Gaussian model is afforded by 
examining the equilibrium (E) as a function of the parameter 13. This is presented in Fig. |4j In the Gaussian noise 
model, the divergence condition condition is obviously satisfied for /3 > /3c, since the right-hand size grows linearly 
with /3. The incipient divergence of (E) for the Gaussian model at /3c = 1.7 is apparent. The binomial model, on the 
other hand, shows no special behavior at the Gaussian critical /3 at which (E) diverges. Rather, (E) exhibits a broad 
minimum at /3 3.9 and then rises toward unity as /3 increases. Even before the divergence, the Gaussian model 
deviates significantly from its binomial counterpoint, since binomial noise is far from Gaussian when N is small. The 
third curve in this figure results from a cutoff version of the Gaussian model, to be discussed later. For smaller noise, 
(equivalently, larger N in the binomial model), however, as depicted in Fig. [s] the difference between the two models 
is small, essentially right up to the divergence, as binomial noise is very well approximated by Gaussian noise, except 
in the tails, which are dynamically irrelevant as long as one is a finite distance below the transition. The onset of the 
deviation for small noise in the Gaussian model is extremely close to the critical /3. For example, (E) crosses zero at 
a distance of the order of 10~^^ from the critical /3 for the parameters of Fig. [sj 



III. SUDDEN/ ADIABATIC (S/A) APPROXIMATION 

Unfortunately, the binomial model does not admit an analytic solution for finite N (which is of course the interesting 
case, since otherwise the effect of the noise is infinitesimal). We can however make use of the natural ordering of 
time scales in the problem to construct a solvable limit. For the circuit to show a significant transient response, 
it is necessary for the time scale of the A dynamics to be much faster than the B dynamics, otherwise the system 
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FIG. 3: (color online). Excerpt of a simulation of the N = 1 binomial model with parameters parallel to those of the Gaussian 
simulation in Fig. [2] fconCo = 1, k^g = 2, a = 1, 7 = 2.5, 5 = 0.4, {3 = 1.7. 
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FIG. 4: (color online). Variation of {E) with /3 for the Gaussian model, Eq. [Tjand the corresponding A'^ = 1 binomial model, 
derived from averaging 10* simulations. The parameters of the Gaussian model (except /3) are as in Fig. [2] for the binomial 
model as in Fig. |3] These results are compared to that of the cutoff model (defined later in the text). 
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FIG. 5: (color online). Variation of {E) with /3 for the Gaussian model, Eq. [Tjand the corresponding binomial model, derived 
from averaging 10* simulations. Here A'^ = 10 in the binomial model, and similarly the Gaussian model has — 4/90, so that 
the Gaussian (E) diverges at /3 = 17.84. The other parameters are as in Fig. |4] These results are compared to that of the 
cutoff model (defined later in the text). 



adapts too rapidly to the changing signal and the transient response is aborted. Furthermore, the time scale of the 
E dynamics, l/(^o + Bq), should be intermediate to those of the A and B fields. If the E dynamics is too fast, 
then the system is limited in any case to respond no quicker than A, and if the E dynamics is too slow, A and B 
have reequilibrated by the time E starts to respond, and again there is no transient response. Also, the time scale of 
the noise dynamics should be intermediate. Too fast noise would just get averaged away, and too slow noise would 
be adapted away. Thus, we are led to consider the limit where the A dynamics is much faster than all the other 
processes, and the B dynamics much slower. This limit is analytically treatable, as we now proceed to show. 

Formally, we define our approximate theory, which we denote the Sudden/Adiabatic (S/A) theory, by taking both 
/3 — >■ and (5 — 0, with a fixed ratio Bp = {(3/6). Since the time scale of the B dynamics is so long, the noise in 
the signal is completely averaged over and we can just set B = Bq = i3p(fconCo/(fconCo + fcoff)), its average value. In 
addition, we take the limit of large a and 7 (with fixed ratio Ap) which guarantees that the activator dynamics is fast 
enough to precisely follow the noise, i.e. A{t) = S{t)Ap. 

We first examine the case = 1. To proceed, we decompose the equation for the probability distribution of E into 
Po.i{E), the joint probability of £' and the input signal S begin or 1, respectively, so that P{E, t) = Po{E, t)+Pi{E, t). 
We immediately derive that, in steady state, 

dP d 

= AjonCoPo - KsPi + ^ ([(So + Ap)E - Ap] Pi) = 

BP B 

= -konCoPo + kosPi + -Q^{BoEPo)=0 (14) 

Adding the two equations and integrating gives 



where Ep = Ap/[Ap + Bq). Substituting this back into Eqs. 14 and defining r+ = kog/[Ap + Bq), r_ = konCo/Ba, we 
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FIG. 6: The probability distribution function, P{E) = P+{E) + P{E) as a function of E for the A*' = 1 binomial model with 
a = 100, 7 = 10, /3 = 0.2, 5 = 0.1, fconCo = k^s = 0.4, compared with the A*' = 1 S/A model, with parameters Ap = 10, -Bo ~ 1- 



can find the normalized probabilities defined on the interval < E < Ep, 

Po = "^tl ^'^ E;(^^+^-^ i^—^ (Ep - (16) 

The total probability P{E) has the interesting behavior of switching from being peaked at the interval center to the 
interval endpoints, as the parameters are varied p3J. For r+ < 1, the probability density diverges at Ep, and for 
r_ < 1, the density diverges at 0. The above expression immediately gives the prediction 



co/(fc 

on Co + fcoff) 

= . ir_ + l ^ ' 



To compare this to the full binomial model, we conducted a simulation with a — 100, 7 = 10 (giving Ap = 10), 
/3 = 0.2, 6 = 0.1, and fconCo = fcoff = 0.4, (giving Bq = 1). Here, the above theory predicts {E)/Eo = .684, where 
the deterministic Eq = fconCo^p/(fconCo^p + (fconCo + fcoff)So). The simulation gave {E) / Eq — .691 and indeed the 
histogram is peaked at the endpoints, as predicted by the analysis. This is seen in Fig. |6] 

In particular, our limiting theory predicts that as a function of Bq, i.e. /?, {E)/Eq starts at a value of unity at 
Bq = 0, which is reasonable since the system is saturated and so (E) is unity independent of the noise. For small Bq, 
(E) / Eq falls with an an initial slope of fcoff/(fconCo(fconCo + fcoff))- As a function of Bq, {E)/Eq reaches a minimum at 
Bq = ^ ApfconCo and then turns back up, approaching unity at large Bq. This qualitative behavior is in accord with 
what we saw in Fig. 1, even though there the parameters are far from fulfilling the separation of scales assumed in 
the analysis. For k^^CQ = fcoff, the value of {E)/Eq at the minimum is 

miEQ)„,,^ = ^(J'"_^'^2r + 2) ' ^ (18) 

This decreases from unity for small r to a value of 1/2 at large r. Thus the larger Ap/ko^CQ, the larger the noise-induced 
relative suppression of {E), since the effective noise amplitude increases as fconCo decreases. When fcoff ^ fconCo, the 
value at the minimum approaches unity, so there is no suppression. On the other hand, when k^s ^ konCQ, the value 
at the minimum approaches fcoff/(fcofi + Ap), which indicates the maximal suppression occurs at Ap ^ fcoS ^ fconCo- 
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A. Moment Equations 

One cannot extend the above approach to compute the an exact closed-form expression for the steady-state distri- 
butions for > 1. However, one can make progress by recognizing that the moment equations take a particularly 
simple form. Consider the steady-state master equation for the case of general N. We have, in an obvious notation, 



= konCoPN-i~2kosPN + ^[[iBo + Ap)E-Ap]P2 



1 ^ 

= -jkonCoPj - (TV - j)koffPj +{N-j + l)fconCoPj-l + (j + l)fcoffP,- + l 
d 

N N 



dE 



Pj] j = 2...7V-l 



= ^NkonCoPo + kosPi + +^{BaEPo) (19) 
Because of the form of these equations, we can get a closed linear system for the moments Zn = J dEPn{E)E: 

= -{N - j)koncozj - jkoffZj + {N-j + l)fco„coZj-i + (j + l)kosZj+i - (^Bo + z^ ~ ^ (20) 

where n„ = / dEPniE) are just the binomial occupation probabilities and z_i = zat+i = 0. This (A^ + 1) x (iV -f 1) 
linear system can immediately be solved for the z's for any given A^. In Fig. [Tj we plot the results obtained by solving 
these equations to determine (i?) = X]^=o 

For large iV, the distribution 11^ becomes highly peaked around its mean, j — N ko^^Ci^ / [ko^c^ -\- fcoff): and so does 
Zj. Thus, to leading order, we can approximate jAp/N by its mean, Aq = Apko^co / {konCo + fcoff), allowing us to solve 
the resulting system, 

z, « zf = Ao/iAo + Bo)n, = EoH, (21) 

Thus, (E) = Eq, and moreover the mean value of E^ conditioned on the value of the input j is in fact independent of j, 
so that adaptation become perfect for large N. To investigate the finite N effects, we expand around Zj, Zj = z r+ Ab- 
using the fact that for all the important modes, j/N — A;onCo/(fconCo + ^off) is small, of order 0(A^~^/^). To solve 
the resultant system, we can approximate it by an ODE for A(y) thought of as a function of the continuous variable 
y = ^/~N{n/N — fconCo/(fconCo + fcoff)- This calculation is presented in Appendix 2, and leads to the result 

AlB^a^ 

" ^" - N{A, + B,nA, + B,+u:) ^^2) 

where w = l/r = /conCo + fcoff- Thus, 

^ " Ar.T(Ao + So)(^o + So+c<.) ^ ^ 

Thus again {E) /E^) starts at 1 for Bq — Q but here falls initially with slope (1 — x)I{Nx{Aq + a;)). Note that this is 
quite different from the B — Q slope for A^ = 1, which was independent of Aq. Again, {E)/Eo has a minimum as a 
function of Bq, here at Bq = Aq^/{1 +uj/Ao), with a value at the minimum of (1 — x)/{Nx{l + + uj/Aq)"^). Thus 
again the suppression is largest for /coff ^ fconCo- Note that for small fconCo, the suppression appears to grow without 
bound, although it is in fact bounded above by unity, indicating that the smaller x, the larger A^ has to be in order 
for the large A^ results to be valid. Nevertheless, the suppression is still strong in this limit. 

IV. GAUSSIAN MODEL - LARGE N LIMIT 

We have seen that our approximate binomial model becomes Gaussian in the large A'' limit. It is interesting to 
check that this agrees with the appropriate limit of the Gaussian model for small noise. For small noise, the Gaussian 
model gives a finite answer, since our criterion is automatically satisfied. Indeed, upon expanding Eq. [7] to linear 
order in (t| and performing the integral we get 



{E) « Eo 



2 SYr(3i-/ - S){-/S + SoP + Soa){aSoST + SoPiT + S^-/T + SYt + jS) 

1 - '^s- 



{aSoS + Sop-f + -/6^){6t + + 6){aSoS + SaPi + -^''5){aSa5T + SoP-fT + 75)(7t + l)5o 



(24) 
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FIG. 7: The scaled suppression of the equilibrium (E), N{1 — (E)) as a function of Bo in the reduced model, for various A'^ = 1, 
2, 5, and 20 along with the asymptotic large-A'' result. The parameters are Aq = 10, fconCo = k^g = 0.4. 



and then taking a and 7 to oo with a/j = Aq/x and taking /? a nd S to zero with P/5 = Bq/x^ and setting Sq — x. 



ag = (J /N, we indeed reproduce our large TV result, Eq. 23 In addition, the general result confirms that the 
suppression is maximized in our distinguished limit 7 3> ^. 

Thus, the Gaussian model is perfectly acceptable in the small noise limit. One can ask if there is way to extend 
it beyond this limit. Clearly just increasing the noise amplitude leads to problems, as we have seen. We have seen 
in Fig. |4] that the problem is not restricted just to noise levels bigger than the critical value. Rather, for this case 
of large noise, the Gaussian answer is not accurate even when we are not close to the critical value. The problem is 
of course the already demonstrated large negative excursions. For small or intermediate noise levels, the problematic 
negative excursions of A and B are actually quite rare. To understand this better, consider the distribution of 
E{t), for some given t. For short times this is well-behaved, but one exceed the time-scale of the E dynamics, the 
distribution develops a power-law tail for large negative E(t)] this can be seen in Fig. [sj If /3 < /3c, the exponent 
of this distribution is greater that 2 in magnitude, and so the first moment is finite. Of course, there is a range for 
which the first moment is finite but the second and higher are already divergent. Since the exact formula shows that 
there is no divergence at finite i, there must be a cutoff in the power-law tail at some extremely large value, which 
however is very difficult to see from the numerics [TSj. For /3 > /3c, on the other hand, the power decreases and the 
first moment diverges as well (subject to the same extremely large cutoff). As the noise level decreases, the above 
picture still hold. The probability of E{t) < 0, however, decreases exponentially with N. Thus, while for /? > /3c, the 
first moment diverges, it becomes exponentially more difficult to see this in a simulation at intermediate noise levels. 
For small noise, it is well-nigh impossible. Thus, simulating the Gaussian model will, for small and intermediate 
noise give perfectly physical answers, which is, however, the incorrect answer for the true ensemble average, which is 
dominated by extremely rare huge events. Of course, if one wishes to rely on simulations, one can simulate directly 
the binomial noise model. 

This line of reasoning leads to an alternate approach for extending the Gaussian model, loosely motivated by the 
successful use of simple cutoffs in regularizing reaction-diffusion equations which arise from the large N limit of Markov 
processes and which overemphasize the growth at small concentrations by missing the essential role of particle number 
discreteness [TBtilSj . (In fact, it has been recently proven that adding a cutoff to the Fisher equation (T^HST] exactly 
reproduces the anomalous front velocity correction for asymptotic large N [22\). Here, we cut off the integral in Eq. 
[7] at the point where either the integrand becomes (unphysically) negative or (unphysically) increasing as u increases. 
Fig. [4] presents results for = 1 showing that this method not only prevents blow-up (by construction) but also 
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FIG. 8: (color online) Gaussian Model: Probability Distribution Function for 1 — E{tf), conditioned on E{tj) < 0, plotted 
in log-log scale, for a — 1.4, 2.2 and 3.0. The tails approach straight lines for large |_E(t/)|, indicating asymptotic power-law 
distributions with approximate exponents of 2.2, 1.8 and 1.6, respectively. A'' — 0.5, t/ — 20 All other parameters are as in 
Fig. 1. For each data set, 10® runs were performed, yielding 146,000, 203,000 and 234,000 data points satisfying E{tf) < 0, 
respectively and the data was binned logarithmically using 25 bins. 

does a decent job in capturing the true answer. It is clear that this method becomes more and more accurate as N 
increases, since for large TV the unphysical neglected integrand is exponentially small in A'^ (even though diverging 
as u — > oo); this can be seen in Fig. [5] . This is reminiscent of the situation we encountered with the Gaussian 
simulations, where the part of the distribution that has diverging mean has exponentially small probability. 

V. DYNAMIC BEHAVIOR 

Up to this point, we have focused exclusively on the steady-state properties of the model. It is also interesting to 
investigate the dynamic behavior of the model, especially since in the deterministic limit adaptation implies that the 
steady-state behavior is independent of the environment, and the only nontrivial behavior is the transient response 
of the system to changes in the input. In our distinguished limit where the A dynamics is fast and the B dynamics is 
slow, we can get a complete picture of the dynamics. Imagine that the input signal suddenly changes to a new value. 
The A field will immediately respond to this change. The B field will only relax slowly to the change. As before, the 
relatively fast noise is averaged over, and so the B dynamics can be taken to be a deterministic exponential relaxation 
to its new equilibrium value Bq, on the slow time scale. Thus, we can consider the changes in B as adiabatic, with 
the system in quasi-equilibrium with the current value of Bq. Thus, the value of E, averaged over the relatively quick 
fluctuations will be that given by our above solution for (E) with Bq = Bo{t) ^ B'q + {B'^ - Bo(ip))e"''(*"*''\ where 
tp is the time of the pulse, and fconCo is given by its new value, with Ap and fcoff unchanged. We see this plotted for 
the case of = 1 in Fig. |9j together with the full simulation, averaged over 10'* runs. We see that this treatment 
captures very well the dynamics on the long 1/6 time scale of the B dynamics. We see that for this strong a noise, the 
transient response of the system is swamped by the shift in the equilibrium value of (E). One can in fact do better 
by solving the time-dependent moment equations, with Ili{t) = x' — [x' — x)e~^^*~*''\ This resolves the initial period 
of the rise in {E{t)), during which the noise equilibrates to its new statistics. This is also shown in the figure. 
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FIG. 9: The ensemble-averaged {E{t)) for the binomial model when cq is suddenly increased by a factor of 2 at = 100. The 
average over 10* runs is denoted by the curve labeled Simulation. The adiabatic approximation where the equilibrium value of 



E is shown with Bo taken to be its ensemble average BO{t) — B'q — {B'q — Bo{tp))e 



-S{t-tp) 



and the other parameters post-pulse 



are taken to be their new values is shown by the trace labeled Adiabatic. The numerical solution of the moment equations 
with Bo{t) taken as above and ni(t) = x' — {x — z)e"''*'~*p' is marked by ODE. The parameters are = 1, a = 100, 13 = 1, 
7 = 10, S = 0.1, and feoff = 0.8. The initial value of feonCo = 0.08, and post-pulse feonCo = 0.16. Thus x = 1/11 and x = 1/6. 
The deterministic equilibrium value of E is 1/2, both pre- and post-pulse. 



VI. SUMMARY 



In summary, we have shown that the Gaussian noise model fails in principle beyond a critical value of /3, the strength 
of the inhibitor production. This renders a full analytic treatment impossible. The cause of this can be traced to 
the interaction of the signaling nonlinearity with the fact that Gaussian noise does not respect the constraint that 
chemical signals must be positive. One can get rid of this effect by linearizing the reaction equation, but this completely 
eliminates the possibility of investigating the extent to which the perfect adaptation seen in the deterministic limit is 
undermined by signal stochasticity. 

Even though it cannot be solved exactly, one can make analytic progress for the most interesting range of parameters, 
namely where the activation is fast and the inhibition slow, as compared to the noise time scale . These can be used 
as a guide to determining the actual deviation of (E) from its signal-independent mean-field value. Also, it is possible 
to devise a cutoff procedure which accurately predicts the results of the binomial model for intermediate and large N 
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Appendix A: Appendix I: Gaussian Noise 

We present here the formulas for $ and AE{u) appearing in the expression for the average E for the Gaussian noise 
model. 

a B 

\ — . , 

7 + l/r 5+ l/r / Vl " 1/ 



1/tJ \j-1/t^ S-1/tJ \ J 



..(^.0% (Al) 



7+l/rV7-l/T ^-l/ry/V 7 



(A2) 



Appendix B: Appendix II: Large TV Limit 

We begin by writing the exact equation satisfied by = Zj — EoTlj: 



= -{N-j)konCoAj-jkosAj+{N-j+l)konCoAj_i+{j+l)kosAj+i~{Bo + 



N 



Aj + ^——A,{l-Eo)Uj (Bl) 



where x = fconCo/ (fconCo + ^off)- To proceed, wc transform this difference equation into an ODE in terms of the variable 
y={j- Nx)/y/N, and writing Aj = N^^ z''^^\y) + N'^/"^ z^'^\y). To leading order in N we get (w = fconCo + kos): 



3)^2770-2 



ye 



the solution of which is 



ApBo 

{Ao + Bo+ w)(Ao + Bo)V2^ 



ye 



(B2) 



(B3) 



This correction refiects the breakdown of perfect adaptation for finite N, but due to its asymmetry, does not lead to 
a correction in {E). To obtain the leading order correction for this quantity, we have to go to next order: 



The solution is then 

z^'\y) = - 



(1) _ ojjl - 2x) (P ^j^) ^ y^{l-2x) _ y^'^ ApBp ^-^ygq^ 



ApBoe-yy^-' 



Apa^ 



dy 



y 



2cr2\/2^ 



3ct2; Aq + Bo 



{Ao+Bo + 2ij) Vf72 ^0 + ^0 



2a; 



(^0 + Bo){Ao + Bo+ uj)V2^ 
The second term does not contribute to {E), interestingly enough, and 

AjBoa^ 



(E) « Eo 



N{Ao + Bo)^{Ao + Bo + ij) 



(B4) 
(B5) 

(B6) 
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